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Abstract 

Perturbation theory is an important tool in the analysis of oscillators and their response to 
external stimuli. It is predicated on the assumption that the perturbations in question are "suffi- 
ciently weak", an assumption that is not always valid when perturbative methods are applied. In 
this paper, we identify a number of concrete dynamical scenarios in which a standard perturba- 
tive technique, based on the infinitesimal phase response curve (PRC), is shown to give different 
predictions than the full model. Shear-induced chaos, i.e., chaotic behavior that results from the 
amplification of small perturbations by underlying shear, is missed entirely by the PRC. We show 
also that the presence of "sticky" phase-space structures tend to cause perturbative techniques to 
overestimate the frequencies and regularity of the oscillations. The phenomena we describe can 
all be observed in a simple 2D neuron model, which we choose for illustration as the PRC is 
widely used in mathematical neuroscience. 

Introduction 

Rhythmic activity is commonplace in biological phenomena: the spontaneous beating of heart cells 
in culture ifTOll . the synchronization of flashing fireflies ll27l . and central pattern generators in ani- 
mal locomotion 131, and calcium oscillations that underlie a plethora of cellular responses (ranging 
from muscle contraction to neurosecretion) [32] are just a few examples (see, e.g., [11] and [38] for 
many more). Mathematical models of biological oscillations often provide useful insights into the 
underlying biological process; for example, they can explain the observed robustness of circadian 
rhythms 1381 and of population cycles 1251 . and can be used to infer plausible structures for central 
pattern generator networks based on locomotion gaits fV2}\. Analyzing models of biological oscilla- 
tions, however, is generally not easy: the mechanisms underlying biological rhythms are varied and 
complex, and this complexity is reflected in the corresponding mathematical models. Indeed, models 
of biological oscillators are often high dimensional, highly nonlinear, and have uncertain parameters, 
all of which make them challenging to study. 

Perturbation theory, long a staple of applied mathematics, provides a practical solution in many 
situations. Mathematically, robust oscillations correspond to attracting limit cycles in phase space. If 
the stimuh involved are not too strong, then one is justified in viewing trajectories of the forced system 
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as perturbations of the original limit cycle. When judiciously applied, such perturbative analyses can 
yield a great deal of insight and useful quantitative predictions. For oscillators, a good example of an 
effective perturbation technique is that of infinitesimal phase response curves (PRCs) and associated 
phase reductions 171. The infinitesimal PRC of an oscillator records the phase change that results 
from an applied perturbation. Given an oscillator, there are many ways to obtain its infinitesimal 
PRC: in addition to analytical perturbation techniques, there are efficient numerical methods for con- 
structing PRCs; the entire PRC itself can even be inferred directly from experimental data. Moreover, 
PRC -based techniques require only tracking just the phase of the oscillator, thereby greatly reducing 
complexity. They are used in many areas of mathematical biology, but especially in computational 
neuroscience, as they yield direct predictions about the modulation of spike timing and frequency 
by external stimuli. Furthermore, they allow one to make predictions for weakly-coupled networks 

nsiiaii. 

However, as is well known, perturbation methods do not always correctly reflect the true dynam- 
ical picture, as they systematically overlook certain aspects of the dynamics. In using perturbation 
theory, one assumes that the perturbation is small, an assumption that is not always valid in applica- 
tions. In this paper, we identify a number of concrete scenarios in which PRCs and phase reductions 
give predictions different from that of the full model. One situation is when the perturbation causes 
the trajectory to leave the basin of attraction of the limit cycle, which can occur even with moderately 
weak forcing. But even without leaving the basin, more subtle effects can lead perturbation theory 
astray, giving — to varying degrees — incorrect predictions. These situations include the presence 
of "sticky" invariant phase-space structures near a limit cycle, which can cause perturbation theory to 
overestimate the regularity and frequency of a stimulated oscillator We will also show that dynam- 
ical shear in a neighborhood of the oscillator can cause it to behave chaotically when forced. These 
scenarios cannot be captured by infinitesimal phase reductions. 

The scenarios described in this paper are relevant for general oscillators, but in view of the pop- 
ularity of PRC -based techniques in neuroscience, we will illustrate our ideas using the Morris-Lecar 
(ML) neuron model HHH. This widely-used model provides a convenient and flexible example be- 
cause of its low dimensionality and rich bifurcation structure. The phenomena of interest are not 
hard to find in certain ML regimes; they do not require stringent tuning of parameters. We point out 
also that while we focus here on single neuron models, our findings remain relevant for oscillators 
operating within networks. 

This paper is organized as follows: In Sect.[Tl we review some relevant mathematical background, 
including brief discussions of phase response theory and "shear-induced chaos", a general mechanism 
for producing chaotic behavior in driven oscillators. Sect. |2] introduces the ML model and some 
relevant ideas from computational neuroscience. In the last three sections, certain regimes of the ML 
model are used to demonstrate how perturbative techniques sometimes do not correctly predict the 
behavior of the full model. Sect. [3] contains an example in which the infinitesimal PRC gives no hint 
of the strange attractor in the full model. Sects. |4]and[5]illustrate how the presence of nearby invariant 
structures can impact neuronal response in ways that cannot be captured by PRCs alone. 

1 Mathematical background 

The general setting for this section is a nonlinear oscillator modeled by an n-dimensional ODE x = 
f{x) with a limit cycle 7 . We assume throughout that 7 is not only attractive as a periodic orbit but 
hyperbolic, i.e., its Floquet multipliers have absolute values < 1. The period of 7 is denoted Per(7). 
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1.1 Phase response curves and phase reductions 

This section contains a brief review of a perturbation technique for oscillators known as (infinitesimal) 
phase response theory. For more details and many applications, see, e.g., ll9l lTniT3l[38l . 

First, we fix a notion of "phase" on the oscillator: Fix a reference point G 7, and declare its 
phase to be ipi^*) = 00 For other point x G 7, the phase 'ip{x) can be defined to be the amount 
of time to go from to x along the cycle 7 . Next we extend this notion to a neighborhood of 7. It 
is a mathematical fact that Tp extends uniquely to a smooth function such that for all solutions x with 
x(0) near 7, ^il){x{t)) = 1 . Thus ijj serves as a kind of "clock" for tracking the passage of time 
along trajectories near 7. 

Our goal is to describe what happens to the phase when a time-dependent perturbation is applied. 
Let 9{t) := ilj{x{t)), where x is a perturbed trajectory. We would like an approximate equation for 
9{t) . Rather than doing this for general perturbations, we consider a perturbed equation of the form 

x = f{x) + I{t)k, (1) 

where x{t) G , I is a scalar signal, and k G M" is a constant vector. (In neuroscience applications, 
for example, one of the variables usually represents the membrane voltage, and this is typically the 
only variable that can be directly affected by external perturbations.) Now let ^ : R ^ R" be a 
periodic solution to the "adjoint equation" 

i{t) = -Df[^{t))^ -m. (2) 

where 7(t) denotes an orbit parametrizing the cycle 7 , and let A(0) := ^{9) ■ k . Under the normal- 
ization condition ^(0) • /(O) = lUit can be shown (see, e.g., |9|) that if x is a solution of Eq. ([T]) and 
/ = 0(e) for a small parameter e, then 6 satisfies 

9 = 1 + A{9)I{t) + 0{e^) . (3) 

Truncating all terms of O(e^) in Eq. Q yields an equation for 6, the phase reduction of Eq. ([Hi. The 
function A is the infinitesimal phase response curve (PRC)I 

For systems that are near a bifurcation, the above procedure can be used to derive analytical 
approximations of the PRC via normal forms Q. In more general situations (where there are few 
practical analytical techniques), one can obtain PRCs through numerical computation, or even directly 
from experimental measurements (see, e.g., [29ll2ll[3B] and references therein). This flexibility and 
accessibility is part of its appeal in mathematical biology, and in neuroscience in particular. 

Stochastic forcing. The basic methodology of infinitesimal phase response curves can be extended to 
systems driven by stochastic forcing. That is, suppose that in addition to a deterministic forcing, we 
add a second white-noise term: 

X = /(x) + {I{t) + /3W) k ; W = white noise. 

In |[24ll . Ly and Ermentrout show that with the above forcing, 9 satisfies the equation 

9 = 1 + A{9)I{t) + y A(e)A'(0) + I3A{9)W + (higher order terms) , (4) 



'in neuroscience, it is customary to define zero phase to be the moment the neuron "spikes.' 
^It is easy to checlc that if this condition holds for t = 0, it holds for all t. 
Some authors refer to A as the phase resetting curve. 



3 



assuming both e and /? are small. 

Eq. dH) can be used to derive a number of quantities of interest. For example, if we view Eq. dH) 
as modeling a neuron that "fires" whenever ^ = , then a result of Ly and Ermentrout states that the 
firing rate resulting from a constant forcing I{t) = e is 



to leading order in /? . In the above, A = A{6) dO. 

Remark on terminology. Other variants of the PRC exist for finite-size perturbations. In this paper, 
the term "PRC" will always mean the infinitesimal phase response curve. 

1.2 When do phase reductions work, and when might they fail? 

Because the PRC is defined in terms of a vector field / and its derivative along a limit cycle 7 (see 
Eq. it can only contain information about the flow in an infinitesimal neighborhood of 7. We 
discuss briefly in this subsection a few scenarios in which the behavior of the flow a finite distance 
away from 7 can have a dramatic effect on the oscillator's response. These ideas are illustrated in 
concrete examples in Sects. |3]-|5] 

(1) Leaving the basin of ^. The simplest possible way for something to go wrong is when the forcing 
carries a phase point outside of the basin of attraction of 7. If denotes the unforced flow, then the 
basin of 7, denoted Basin(7), is defined to be the set of all phase points x such that <l*t(2;) — )• 7 as 
t ^ 00. Once a perturbation, say in the form of a "kick", causes a trajectory to leave Basin(7), what 
it does may depend on dynamical structures far away from 7. For example, if the system is bistable, 
or multi-stable, i.e., it has more than one attracting set, which can be in the form of a stationary point, 
a limit cycle, or something more complicated, then an "escaped" trajectory can end up near one of 
these structures, and in time, it may - or may not - get kicked back into Basin(7). Needless to say, 
the behavior of such a trajectory bears little resemblance to that predicted by the PRC. This scenario 
must be taken into consideration when the forcing is strong relative to the distance of 9(Basin(7)) to 
7. (Here (9(Basin(7)) is the boundary of Basin(7).) 

(2) Invariant structures and "trapping" . Even without venturing outside of Basin(7), a perturbed 
trajectory that comes near (9(Basin(7)) can be nontrivially affected by certain dynamical structures 
in (9(Basin(7)). These structures may seem innocuous - they are non-attracting - but as we will see, 
they can seriously impact the surrounding dynamics. Consider, for example, a saddle fixed point of 
the unperturbed flow. An orbit that comes near it will, under the unperturbed flow, remain near it 
for some duration of time (depending on the ratios of the eigenvalues of the linearized flow at that 
point). Since an invariant structure is typically not truly invariant for the perturbed flow, a perturbed 
trajectory that gets near it will likely escape eventually and return to 7. However, the escape time can 
be quite long, and this effect cannot be captured by phase reductions. The tendency to remain near 
invariant structures can be mitigated by the forcing when the forcing acts to push trajectories away; 
by the same token, it can also be magnified if the forcing "conspires" to keep trajectories in a region. 
Indeed, there is no reason why a forcing cannot - by itself - create trapping regions within the basin 
of 7 if the attraction to 7 is weak compared to the forcing. 
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SHEAR 




Figure 1. The stxetch-and-fold action of a kick followed by relaxation in the presence of shear. 

(3) Shear-induced chaos. This is yet another dynamical phenomenon that cannot be captured by 
(infinitesimal) PRCs. This phenomenon is illustrated in Fig. [T] In each of the two pictures, 7 is 
represented by the horizontal line. Shear refers to the differential in the horizontal component of the 
velocity as one moves vertically up in the phase space; here points above 7 move around 7 faster than 
points below. Suppose that an impulsive perturbation, or "kick," applied in such a situation produces 
a "bump" in 7. As the flow relaxes, this "bump" is attracted back to the original limit cycle. As it 
evolves, it is folded and stretched by the flow if sufficient shear is present, as one can visualize in 
Fig. [T] Since stretching and folding of phase space is associated with complex dynamical behavior 
such as horseshoes and strange attractors |[T4l . this picture suggests that perturbing a limit cycle 
with strong shear can lead to chaotic dynamics. That this can indeed happen has been established 
(rigorously) in recent developments in dynamical systems theory. 

To connect paragraphs (2) and (3), we mention that invariant structures in 5(Basin(7)) can be 
a contributing factor to shear, but shear can also arise for many other reasons. Since the results on 
shear-induced chaos alluded to above are not widely known in the mathematical biology community, 
we will provide a more detailed review in the next subsection. 

1.3 Shear-induced chaos and related phenomena 

In this section, we review some of the geometric ideas put forth in ll35l[36l (and also ll34l[37]| and |[30]| ). 
The exposition here roughly follows |19|, which contains a more thorough (but still non-technical) 
discussion. We focus here on periodically -kicked oscillators, since in this setting the various dynam- 
ical mechanisms are most transparent. 

A periodically-kicked oscillator is a system of the form 



where {A, T) € M x M+ are parameters and H : R" R" is a given smooth function. We assume as 
before that x = f{x) has an attracting hyperbolic limit cycle 7 . Eq. (O thus models an oscillator that 
is given a sharp "kick" every T units of time. We interpret the kicks as follows: whenever t = nT, 
we apply a mapping k (defined by H) to the system; between kicks, the system follows the flow 
generated by x = f{x). When kicks are applied repeatedly, the dynamics of Eq. Q can be 
captured by iterating the time-T map Ft = ° k,- If there is a neighborhood U of ^ such that 
k{U) C Basin(7), and T is long enough that points in k{L() return to U, i.e., Ft{U) C U, then 
r = r\n>oFf{U) is an attractor for the periodically kicked system Ft. One can view F = T{k, T) as 
what becomes of the limit cycle 7 when the oscillator is periodically kicked. 



00 
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The structure of T and the associated dynamics depends strongly on the kick parameters A and 
T, as well as on the relation between the kick map k and the flow near 7 . When A is small, we 
generally expect F to be a slightly perturbed version of 7 , because (as is well known) hyperbolic 
limit cycles are robust under small perturbations |[T4ll . In this case, F is known as an invariant circle, 
and the restriction of Ft to F is equivalent to a diffeomorphism on . Circle diffeomorphisms 
are well known to exhibit essentially two distinct types of behavior: quasiperiodic motion, in which 
the mapping is equivalent to rotation by an irrational angle, and gradient-like behavior characterized 
by sinks and sources on the invariant circle. In terms of the kicked oscillator dynamics, the former 
corresponds to the driven oscillator drifting in and out of phase relative to the kicks, while the latter 
corresponds to stable phase-locking. 

The preceding discussion suggests that when kicks are weak, we should expect fairly regular 
behavior. To obtain more complicated behavior, it is necessary to "break" the invariant circle. The 
main idea is best illustrated in the following linear shear model, a version of which was first studied 
in BOl : 



where {9, y) € x M are coordinates in the phase space. A, a, and A are constants, and : — )• R 
is a non-constant smooth function. When A = 0, the unforced system has a limit cycle 7 = 5^ x {0}. 
The following result, due to Wang and Young, shows that Eq. ^ indeed exhibits chaotic behavior 
under the right conditions. (There is an obvious analog in n-dimensions 1,36,1 .) 

Theorem 1.1. [.36,1 Consider the system in Eq. (fT]). If the quantity 



is sufficiently large (how large depends on the forcing function H), then there is a positive measure 
set A C M"*" such that for all T A,T is a "strange attractor" of Ft ■ 

It is important that H be non-constant, as H{6) is what creates the bumps in Fig.[T] The geometric 
meaning of the term involving cr, the shear, is as depicted in Fig.[T] It is easy to see why ^ • A is key 
to production of chaos by fixing two of these quantities and varying the third: the larger cr, the larger 
the fold; the same is true for larger kick size A. Notice also that weaker limit cycles are more prone 
to shear-induced chaos: the closer A is to 0, the slower ^(7) returns to 7, and the longer the shear acts 
on it, assuming T is large enough. 

The term "strange attractor" in Theorem 11. 11 is used as short-hand for an attractor with an SRB 
measure^ which roughly speaking, implies that the trajectory is unstable, or has a positive Lyapunov 
exponent, starting from Lebesgue-almost every initial condition in the basin of the attractor (or at 
least in a positive measure set). We say such a system has a "strange attractor" because it exhibits 
sustained, observable chaos, i.e., chaotic behavior that is sustained in time, and observable for large 
sets of initial conditions. This is a considerably stronger form of chaos than the presence of horseshoes 
alone (see e.g. ^141 for a discussion of horseshoes). In the latter, it is possible for almost all orbits 
to head toward stable equilibria, resulting in negative Lyapunov exponents; this scenario, in which 
horseshoes coexist with sinks, is known as transient chaos)^ 

''For more information, see |6l l39| . 

'Note that Theorem I 1 . 1 [ asserts the existence of "strange attractors" only for a positive measure set of T, not for all large 
T . Indeed, there exist arbitrarily large T in the complement of A for which Ft exhibits only transient chaos. 
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Figure 2. Geometry of folding in relation to the M^'*'* -foliation. Shown are the kicked image of a segment 
70 and two of its subsequent images under , p = Per(7). 

Shear-induced chaos and the geometry of strong stable manifolds or isochrons 

We now return to the general setting of Eq. ^ and seek to understand what plays the role of the 
shear (as a is no longer defined). Let 7 and $4 be as at the beginning of Sect. 11.31 Crucial to this 
understanding is the following dynamical structure of the unperturbed ilow '■ For x G 7, define the 
strong stable manifold or isochron through x to be 

W{x) = {y : \^t{y) - ^t{x)\ ^ as t ^ 00} . 

With 7 assumed to be hyperbolic, it is known that (see, e.g., |[T3l ) 

(i) W'^^{x) intersects 7 transversally at exactly one point, namely x , and these manifolds are 
invariant in the sense that ^tCW^^^x)) = W^^i^tix)) . 

(ii) {W^^{x),x € 7} partitions Basin(7) into codimension- 1 submanifolds. 

We now examine the action of the kick map k in relation to the W^*''-foliation. Fig.|2]is analogous 
to Fig. [TJ it shows the image of a segment 70 of 7 under Ft = ° ^- For illustration purposes, we 
assume 70 is kicked upward with its end points held fixed, and assume T = np for some n E Z+ 
(otherwise the picture is shifted to another part of 7 but is qualitatively similar). Since $np leaves 
each l^^^-manifold invariant, we may imagine that during relaxation, the flow "slides" each point 
of the curve ^(70) back toward 7 along VF^'^-manifolds. In the situation depicted, the effect of the 
folding is evident. 

Fig. IDgives considerable insight into what types of kicks are conducive to the formation of strange 
attractors. Kicks along VK'^^-manifolds or in directions roughly parallel to the VF^'^-manifolds will not 
produce strange attractors, nor will kicks that essentially carry one VF^^-manifold to another. What 
causes the stretching and folding is the variation in how far points x G 7 are moved by k as measured 
in the direction transverse to the VF***-manifolds. In the simple model of Eq. dV]), because of the 
linearity of the unforced equation, VF'^^-manifolds are straight lines with slope — X/a in {6,y)- 
coordinates. Variation in the sense above is created by any non-constant H; the larger the ratio jA, 
the greater this variation. 

The notion of "phases" in Sect. 11.11 is defined precisely by the partition of neighborhoods of 7 
into VF^*-manifolds or isochrons (two names used by different communities for the same object), 
i.e., we view x E Basin(7) as having the same phase as y G 7 if x G W^''{y). The ideas in the 
last paragraph are in the same spirit as the phase transition curves introduced by Winfree ||38l[m i9l. 
They were discovered independently in the rigorous work of Wang and Young, who proved, under 
suitable geometric conditions on phase variations, the existence of strange attractors having many 
of the properties commonly associated with chaos |[34ll37ll . These ideas have since been applied to 
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various situations; rigorous results include Il35l[36ll23l [301 1331 and fSl, and numerical results indicate 
tiie occurrence of siiear-induced ciiaos in broader dynamical settings |[T8ll20l[T7l . 



1.4 Summary and comparisons 

Sects. [TTl and l 1 Jl outlined two seemingly distinct approaches to the analysis of phase response. These 
approaches are in fact closely related: In Sect. ll.ll perturbations are assumed to be small, which geo- 
metrically means one can approximate W^'^{y) by its linearization around 7 . The procedure of phase 
reduction then amounts to viewing the perturbation as a sequence of infinitesimal kicks, projecting the 
kicked trajectory back to 7 along VF*'^ -leaves following each kick. The analysis sketched in Sect. ll.3l 
is a more global version of the same idea: here the perturbed orbit is allowed to wander farther away 
from 7, and one takes into consideration the geometry (or curvature) of the l^/'^^-manifolds in relation 
to the kick in assessing its impact. 

A notable difference between full-model analyses and phase reductions is that the latter rule 
out a priori any possibility of chaotic behavior. As explained in Sect. II. 3[ in the full model, fairly 
innocuous kicks applied "the right way" can lead to positive Lyapunov exponents for large sets of 
initial conditions. This cannot be captured by the infinitesimal PRC, for flows in one spatial dimension 
are never chaotic. 

With regard to practicalities, full model analyses are, needless to say, more costly. While our aim 
here has been to raise awareness of the issues discussed in Sects. [L2l and ll.3[ it is not our intention to 
advocate that one necessarily starts by computing VF*''-manifolds (even though they are computable 
in many situations). In many cases, by far the most direct way to get a quick idea of whether shear- 
induced chaos is present is to look at Fy-images of 7 (recall that Ft = o k is the composite 
map obtained by first kicking then following the unperturbed dynamics for time T) and see if folds 
develop. See Section 3. 

2 A neuron model 

In Sects. [3]-[5] a few concrete scenarios will be presented to illustrate how and why the infinitesimal 
PRC may give incorrect predictions. The model we use is taken from neuroscience; it is the Morris- 
Lecar (ML) model of neuron dynamics. A brief introduction of the ML model and some relevant 
information is given below to provide context for readers not familiar with the subject. We remark 
that perturbation methods and the infinitesimal PRC in particular are widely used in neuroscience, 
both in the study of single neuron dynamics (see e.g. 191) and in the analysis of neuronal networks 
modeled as systems of coupled phase oscillators (see e.g. |[T6l ). 

The Morris-Lecar (ML) model 

The dominant mode of communication between neurons is via the generation and transmission of ac- 
tion potentials, or "spikes" [4J. Neurons accomplish this through the coordinated activity of voltage- 
sensitive ion channels in the cell membrane, which open and close in specific ways in response to 
changes in membrane voltage. The ML model is a simple model of this spike-generation process for 
a single neuron. It has the form 



(8) 
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The variable v is the membrane voltage; the first equation expresses Kirchoff 's current law across the 
cell membrane, with I{t) representing a stimulus in the form of an incoming current. The variable 
w is a. gating variable: it describes the fraction of membrane ion channels that are open at any 
time. Real-life neurons typically have multiple kinds of ion channels; in more realistic models like 
Hodgkin-Huxley, these are tracked by separate gating variables. The ML model is simplified in that 
there is just one effective gating variable. The forms and meanings of the auxiliary functions u^oo > 
Tw , m-oo and other parameters in ([8]) are given in the Appendix. For more information we refer the 
reader to IH. 

The ML model has a very rich bifurcation structure. Roughly speaking, by varying a constant 
current I{t) = Iq , one observes, in different parameter regions, dynamical regimes corresponding 
to sinks, limit cycles, and Hopf, saddle-node and homoclinic bifurcations, as well as combinations of 
the above. These scenarios together with their neuroscience interpretations are discussed in detail in 
|[9|| . We will use two of these scenarios for illustration; relevant features of these regimes are reviewed 
as needed in the sections to follow. 

Oscillators in neuronal dynamics 

When a sufficiently large DC current is injected into a neuron, either artificially or through the action 
of neurotransmitters, a typical neuron will begin to spike regularly. In such situations, one can view 
the neuron as an oscillator. For example, in Eq. if we apply a constant driving current I{t) = Iq 
and slowly increase Iq from to a large value, the ML neuron will switch from quiescent to spiking 
at regular intervals, the latter corresponding to the emergence of a limit cycle. More generally, if a 
neuron is operating in a "mean-driven" regime, in which the stimulus it receives consists of a large 
DC component plus a (weaker) fluctuating AC component, one can view the AC component of the 
stimulus as a perturbation of the oscillator ||9l. 

Relevant properties 

The crudest statistic associated with a spiking neuron it is firing rate, which translates into the fre- 
quency of the neural oscillator. Sometimes one is also interested in more refined properties of spike 
trains and even the precise timing of spikes. 

A neuron or a network of neurons receiving a stimulus is said to be reliable if its response does not 
vary significantly upon repeated presentations of the same stimulus. Reliability is of interest because 
it constrains a neuron's (or network's) ability to encode information via temporal patterns of spikes. 
Mathematically, a stimulus-driven system can be viewed as a non-autonomous dynamical system of 
the form x = f{x, I{t)) , where I{t) represents the stimulus. The question of spike-time reliability, 
then, boils down to the following: Given a specific signal (/(t) : t G [0, oo)) , does the response x{t) 
depend (modulo transients) in an essential way on x(0), the condition of the system at the onset of 
the stimulus? If the answer is negative, the system is reliable. Otherwise, it is unreliable. 

The relevant dynamical quantity here is Amax. the largest Lyapunov exponent of the system : If 
Amax < 0, then the system is reliable, whereas Amax > leads to unreliability. Heuristically, this 
is because A^ax < leads to phase space contraction, so that the effects of initial conditions are 
quickly forgotten, whereas A^ax > leads the system to amplify small differences in initial states. 
The reasoning can be made more precise via the theory of random dynamical systems and random 
attractors; see |l20l|2T]|22]| for details. 
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(a) At the moment of a homoclinic bifurcation. Iq 35 (b) The "homoclinic regime" at Iq = 39.5 

Figure 3. A homoclinic bifurcation and the "homoclinic regime." In this parameter regime, the system 
has 3 fixed points (shown as black dots), the middle one being a saddle for a broad range of Ig . Panel (a) 
shows a homoclinic loop anchored at the saddle, with a source on the right and a sink on the left. As Iq 
increases, the homoclinic loop breaks apart and a limit cycle emerges. Panel (b) shows the phase portrait at 
the parameter regime we use, which is well past the homoclinic bifurcation. In addition to the homoclinic 
bifurcation shown in (a), a (subcritical) Hopf bifurcation has occurred, leading the fixed point on the right 
to become a sink surrounded by an unstable periodic orbit / cycle (dashed curve). 

3 Chaotic response to periodic kicking 

In this section, we show numerically that shear-induced chaos occurs in the "homoclinic regime" of 
the ML model (see below), leading to a lack of reliability. As noted in Sect. 11.41 such a possibility is 
ruled out a priori by the infinitesimal PRC. 

3.1 Geometry of the "homoclinic regime" 

We view Eq. ^ with I{t) = Iq for some fixed /q as the unperturbed system, and apply to it a forcing 
in the ?7-variable (forcing the system in w has no physical meaning). 

At Iq = /crit ~ 35, for suitable choices of parameters (details of which are given in the Ap- 
pendix), Eq. dl]) has a homoclinic loop, i.e. there is a saddle fixed point p one branch of whose stable 
and unstable manifolds coincide (see Fig. Oa)). To the left of p lies a sink, which attracts the left 
branch of the unstable manifold; and there is a source inside the loop. 

For Iq > /crit> the homoclinic loop is broken, with the unstable manifold "inside" the stable 
manifold. Fig.[3tb) shows the phase flow at Iq = 39.5. The unstable manifold wraps around a newly 
emerged limit cycle, which will be our 7. Notice also the other dynamical structures: the saddle p, its 
stable and unstable manifolds and the sink to the left, plus a new sink and unstable cycle (indicated 
by the dashed loop) that emerged from the original source via a Hopf bifurcation (the latter occurs 
around Iq « 36.3). In the rest of this section, Iq will be taken around 39.5, and the unforced flow is 
qualitatively similar to that in Fig. Eh). 

We consider periodic kicks applied to such a regime, i.e., in Eq. (HJ we take as input I{t) = 
Iq + AJ2k^i't ~ 1 where \A\ is the kick amplitude and T the kick period. Geometrically, this 
kick corresponds to shifting simultaneously all phase points by A (which can be positive or negative) 
in the horizontal direction. 
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If \A\ is sufficiently large, some very obvious things can happen: For example, a large kick to 
the left can easily drive points on 7 to the left of the stable manifold of the saddle; these points will 
then head for the sink on the left and possibly never return. Another possibility is to drive points 
into the basin of attraction of the sink inside the unstable cycle. We will show in the next subsection 
that something more subtle can happen even with small to medium kicks which do not drive points 
outside of the basin of 7. 

3.2 Shear-induced chaos 

We will proceed in three stages: First, we discuss - on a theoretical level - the dynamical mechanism 
responsible for shear production in the regimes of interest. Then we perform some relatively quick, 
exploratory numerical simulations to check that this shear is sufficiently strong, and to identify suit- 
able parameters. Finally, we confirm the presence of shear-induced chaos via careful computations 
of Lyapunov exponents. 

What is the mechanism that produces shear in the present setting? 

Take A < 0, so that the kick k moves the limit cycle 7 to the left of itself. This pushes about half of 7 
inside the unstable cycle, and half of it outside and to the left. Now suppose we go through a period of 
relaxation, i.e., we apply the ML flow ^t, no kicks. As points come down the left side of 7, those that 
are closer to the stable manifold of the saddle p are likely to follow it longer; consequently they will 
come closer to p. (We assume the kick is small enough that no point gets kicked to the other side of 
the stable manifold.) It is easy to see that the closer an orbit comes to a saddle fixed point, the longer 
it remains in its vicinity - certainly longer than orbits that are, for example, inside the unstable cycle. 
The differential in "time spent near p", if strong enough, may lead to a fold in the o K)-image of 
7- 

The argument above can be made rigorous. Similar ideas have been used in rigorous work lfTl[33l. 
But the reasoning is qualitative: while it shows that some amount of shear is present, it does not tell 
us whether it is sufficient to cause chaos. 

Simulations to detect shear-induced chaos 

As noted in Sect. 11.41 by far the most direct way to detect the presence of shear-induced chaos is to 
plot the images of ^(7) under the unperturbed flow, and to see if folds develop in time. A sample 
"movie" is shown in Fig.lH As one can see, by t = 10, a "tail" has developed: some points in this tail 
are evidently stuck near the saddle, while some other points, evidently <I>t-images of points kicked 
inside the unstable cycle, have remained inside, and at t = 10 they are beginning to gain on points in 
the tail in terms of their angular- position (or phase) around 7. At t = 15, these points have overtaken 
those in the tail, and the difference is further exaggerated in the last two frames. One would conclude 
that for the parameters shown, the system very likely has shear-induced chaos. 

Because this step can be done quickly for the ML model, we use it to locate parameters with the 
desired properties. The kick size used in Fig.|4]is A = —2, which is quite reasonable physiologically: 
it takes at least 10-15 kicks of this size to push a neuron over threshold. Fig.|4]also tells us that it takes 
on the order of 15 units of time for the fold to begin to form, so that for a periodically-kicked system 
to produce chaos, the kick period should probably be upwards of 20 units of time. (Kicks delivered 
too frequently may also drive points to the left of the stable manifold of p; once that happens, it will 
end up near the left sink.) 

Fig. [5] shows the strange attractor that results from a periodically -kicked regime. 
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Figure 4. Shear-induced folding in the Morris-Lecar model. Here, the stable cycle is given one kick with 
A = —2, then allowed to relax back to the cycle. The limit cycle 7 is shown in gray; also shown are the 3 
fixed points (dots), the stable and unstable manifolds of the saddle, and the unstable cycle (dashed curve). 
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Lyapunov exponents 

Let us now explore the chaotic behavior more systematically. Fig. Oa) shows the largest Lyapunov 
exponent A^ax as a function of the kick period T, for a variety of kick amplitudes. More precisely, 
given a phase point {vq,wq) and a tangent vector ryo. we let rj{t) be the solution of the variational 
equation for the kicked system, and define 

Amax(wo,?^o;%) = lim ylog|??(t)| . 

It is a standard fact that for each {vq,wq), Amax(^^o, ^^o; %) is equal to the largest Lyapunov exponent 
except for r/o in a lower dimensional subspace (see, e.g., In many (but not all) of the cases 
considered, this number is independent of the choice of {vq^wq). 

For small A, the exponents are predominantly negative, just as we would expect: the dynami- 
cal landscape is characterized by sinks and saddles. As A increases, the tendency to form positive 
exponents becomes greater, so that for A = —2 and sufficiently large T, most exponents sampled 
are positive, confirming the strong form of (sustained and observable) chaos discussed in Sect. 11.31 
Note, however, that this is only a tendency: the fluctuations seen in the Lyapunov exponents as we 
vary T are likely real, and reflect the competition between the "horseshoes and sinks" and "strange 
attractors" scenarios. 

As explained in Sect. |2l Lyapunov exponents are useful as probes of neuronal reliability. This is 
illustrated in Fig. Ob). The left panel there shows the response of a system with Amax < 0: clearly, 
after some transients, the response of the system is the same across all the trials. In the right panel, 
we have Amax > 0, and the resulting spike times show significant variability across trials. Moreover, 
this variability will persist in time due to the presence of sustained chaos. 

We finish with the following two remarks: (1) While we have focused on the periodically kicked 
case, similar-sized kicks that arrive at random times that are, on average, sufficiently far apart will 
lead to time-dependent strange attractors for much the same reasons as we have discussed. We will 
not pursue this here, but refer the reader to lITSl l20l . (2) Notice that the saddle p, which is the root 
cause of the chaos in the kicked system in the sense that it is responsible for scrambling the order (or 
phases) of the points on ^(7), does not lie that close to the limit cycle. Not only would phase reduction 
methods miss this effect, but the PRC will offer no hint at all that any breakdown has occurred. 

4 Firing rates and interspike intervals 

We demonstrate in this section that, as asserted in item (2) of Sect. ll.2[ "sticky" sets in the boundary 
of Basin(7) can lead to discrepancies between the dynamical behavior of the full model and that 
deduced from the infinitesimal PRC. Specifically, we will present numerical data to show in the 
example considered ihat firing rates for the full model are significantly lower and interspike intervals 
(ISIs) are longer and more spread out than those predicted by the PRC. 

For the unperturbed flow, we continue to use the example from the last section, with the same 
parameters. Recall that this system is in the "homoclinic regime" of ML, with a limit cycle we call 7. 
Shown in dashed lines in Fig.|5]are part of 9(Basin(7)), the boundary of the basin of attraction of 7: 
one component is the stable manifold of p; the other is a repelling periodic orbit. As discussed 
in Sect. 13.21 all orbits close enough to will follow it to p and remain there for some time before 
moving away, {p} is a candidate "sticky set" in this example. The repelling periodic orbit is in some 
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(b) Spike times over repeated trials 

Figure 6. Lyapunov exponents and reliability of the ML system. In (a), we plot the Lyapunov exponents. 
In (b), we plot spike times generated by reliable (left) and unreliable (right) systems over repeated trials, 
with random initial conditions sampled from a neighborhood of 7. The exponents in (a) are computed as 
follows: For each choice of kick amplitude and kick period, 6 random initial conditions are used to estimate 
exponents via long-time simulations. The max and min are treated as outliers and discarded; the remaining 
are used to form error bars. The dot marks the median. 
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Figure 7. Firing rate statistics of the stochastically-forced ML model. In (a), we compare of firing rates for 
the full ML system and its phase reduction driven by white noise. Black dots: Full ML. Red squares: Em- 
pirical phase reduction firing rate. Solid blue curve: Perturbative formula for PRC firing rate (see Eq. Q). 
Panel (b) shows a long sample path for /3 = 0.2 . 

sense also a sticky set, especially if the expansion away from the orbit is weak (it is, since it has just 
emerged from a Hopf bifurcation). Indeed orbits that come close to it may cross over to the other side 
and get pushed toward the sink. 

Instead of periodic forcing, here we drive the oscillator with white noise in addition to the steady 
current Iq, i.e., we let I{t) = lo + P-Cm-wE and consider the random dynamical system that 
results from looking at one realization of at a time. For the phase reduced model, it is natural to 
regard the system as producing one spike as it passes through a marked point on the cycle. For the 
full ML model, it is necessary to fix an artificial definition of what it means for the system to "spike": 
we define a reset to be when the voltage falls below -10, and say the neuron "spikes" when (after a 
reset) the voltage rises above -1-12.5. See Fig.[5]for how this corresponds to the location of the limit 
cycle. 

Comparison of firing rates. The results are shown in Fig. [V^a). Plotted are firing rates as a function 
of the drive amplitude /3 for three different systems: (i) the full ML system; (ii) the firing rate of the 
phase reduction, computed empirically; and (iii) the firing rate of the phase reduction as computed by 
the perturbative formula ^ of Ly and Ermentrout. We plot both (ii) and (iii) because Eq. ^ is itself 
an approximate result based on the phase reduction (0]), valid only for small /3. As one would expect, 
for smaller forcing amplitudes (< 0.1), all three agree. As (3 increases, the perturbative formula tracks 
the empirical firing rate of the phase reduction fairly well, but neither of these PRC -based predictions 
capture the dramatic drop in firing rate of the full system occurring around /? = 0.2. Simulations (an 
example of which is shown in Fig. Hb)) show that the precipitous drop in firing rate in this case has 
to do with trajectories spending more time near the saddle p. 

Comparison of ISI distributions. Fig. [8]shows the numerically computed ISI distributions for the full 
ML model and for the phase reduction (01)0 For small noise, we see the full ML system and phase 
reduction agree fairly well: the ISI distribution is concentrated around the period of the cycle (about 
25.2), and the shape is roughly gaussian. As noise amplitude increases, we see in the full ML system 

*The constant Cm is the membrane capacitance. With this scaling, j3 has the dimension voltage/ \/time ; this makes 
its magnitude relative to the membrane voltage easier to assess. 

^We note that [24) also contains an explicit small-/? expansion of the ISI distribution. We do not use it in Fig.[8]because 
it is not applicable except (possibly) for /3 = 0.05 . 
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Figure 8. interspike interval (ISl) distributions for various /3. Solid black lines: ISI histograms for the full 
ML system. Dashed blue lines: ISl histograms for the phase reduction (|4|. 

(i) a broadening of the ISI distribution, and (ii) an overall shift toward larger ISIs. Neither of these 
effects are captured by the PRC, as they involve structures a finite distance away from 7 . 

Finally, to compare the tails of the ISI distributions, here are the fraction of ISIs greater than 
2 X (cycle period): 



/3 


Full ML 


PRC 


0.25 


0.0225 


0.0003 


0.5 


0.127 


0.0032 


1.0 


0.322 


0.0013 



Notice that the PRC systematically under-predicts the probability of a long interspike interval, con- 
sistent with the "trapping" effect of nearby invariant structures. 

5 Altered spiking patterns in bistable systems 

In this last section, we present a simple example in which the perturbation takes an orbit outside of 
Basin(7), the scenario described in item (1) of Sect. ll.2l It is no surprise that PRC predictions would 
break down here. Perhaps the lesson to take away from this example is that even a weak drive can 
lead to spike patterns that are nontrivially different. 

For the unperturbed system, we consider here parameters of ML that put it in a "Hopf regime", 
with Jo ~ 90. Precise parameters used are given in the Appendix, and some relevant dynamical 
structures are shown in Fig. |9] In this regime, the system has a limit cycle, which we call 7. Each 
time an orbit passes near the right-most point of 7, we think of it as producing a "spike". Notice that 
this system is bistable: inside 7 there is a sink, the basin of attraction of which is separated from the 
basin of 7 by a repelling periodic orbit (from a Hopf bifurcation). This repelling periodic orbit is very 
close to 7, but since it is a positive distance away, it will not show up in PRC considerations. 

We now consider perturbations to the ML model above. As in Sect. HI we again use a white-noise 
perturbation, i.e., I{t) = Iq + /3 • Cm ■ W . Fig. [TO] shows the resulting voltage traces of the full 
model and PRC predictions (computed at Iq = 90). The top panel shows the full ML simulation with 
/3 = 0.2 (a relatively weak forcing), starting the trajectory on 7 . Since the basin of the sink, i.e., the 
unstable periodic orbit shown in Fig. |9l is so close to the limit cycle 7, a trajectory following 7 can 
easily get pushed into the basin, and then attracted to the sink. This can happen even with very small 
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Figure 9. ML phase portrait in the Hopf regime. Shown are a limit cycle surrounding an unstable cycle, 
which in turn surrounds a sink. Here, /o — 90, corresponding to an unforced period of about 102 ms. 
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Figure 10. Voltage traces for the ML system in Hopf regime. The top two panels show simulation using 
the full ML model with the indicated level of stochastic forcing. The bottom panel shows the PRC using 
/3 = 0.4. Here, Iq = 90. 
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/3. The sink itself, however, is a little farther from the boundary of its basin, so it is more difficult for 
a trajectory near the sink to escape under weak /3. This is why for weak noise like /3 = 0.2, it is easy 
to observe a transition from 7 to the sink, but not easy to see transition in the reverse direction. 

For (3 a little larger, e.g., /? = 0.4, the trajectory jumps back and forth more readily: it alternately 
follows (roughly) the limit cycle 7 and stays near the sink, switching between the "spiking" and 
"quiescent" modes at somewhat random times. Not surprisingly, in the phase-reduced model (bottom 
panel), these perturbations do not have a significant impact, and the spiking ranges from completely 
regular to slightly irregular in the second case due to the effect of the noise. PRC voltage traces for 
/3 = 0.2 (not shown) are quite similar. 

Finally, we note that in this bistable situation, a neuron can exhibit substantial sub-threshold 
activity. The PRC underestimates the extent of such activity. 

Conclusions 

This paper compares two ways of evaluating an oscillator's response to perturbations: phase reduc- 
tions versus analysis of the full model. We have found that the infinitesimal PRC, which has the 
virtue of being both straightforward to execute and reducing model dimensionality, produces regular 
oscillatory behavior even when the full model does not. Specifically, we presented examples from the 
Morris-Lecar neuron model to show that 

(i) Periodic kicking of the ML system can lead to unreliable response in the full model via the 
mechanism of shear-induced chaos, contrary to PRC predictions. 

(ii) When stochastically driven, stickiness of nearby invariant structures can lead to lower firing 
rates and longer ISIs compared to PRC predictions. 

(iii) The forcing need not be strong to bring about serious discrepancies in firing patterns between 
full and phase-reduced models. 

Moreover, in all the situations examined, the phase reduction itself offers no hint that any breakdown 
has occurred. 

In terms of neuronal response, our results have the following interpretation: Under certain condi- 
tions, PRCs may overestimate spike-time reliability and firing rates; they may also underestimate the 
mean and variance in interspike intervals, and have a tendency to downplay sub-threshold activity. 
Caution needs to be exercised when interpreting results that come from phase reduction arguments, 
especially for systems near bifurcation points. 

While we have focused on the ML model for illustration, the geometric ideas we used are quite 
general, and we expect similar phenomena to occur in a variety of settings that involve rhythms and 
oscillations. 

Finally, given that pure phase reductions cannot always capture the behavior of high dimensional 
nonlinear oscillators, are there alternative methods that give better characterizations? One classical 
technique, which has seen relatively little application in biological modeling, is that of moving-frame 
coordinates for the analysis of periodic orbits. This is nicely described in Hale [15] and recently advo- 
cated in a neuroscience context by Medvedev [,26.1 . In essence this offers a coordinate transformation 
to a phase-amplitude system that allows one to track the evolution of distance from the cycle as well 
as phase on the cycle. We are currently developing this approach and using it to better understand the 
three points listed above. This, and a framework for understanding the dynamics of weakly coupled 
phase-amplitude models, will be presented elsewhere. 
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Appendix. Morris-Lecar model details 

Here we briefly summarize the details of the ML model used in this paper The interested reader is 
referred to 181191 for more details. 
Recall the ML equations 

CmV = I{t) - g\eak- {v -viei,k) - 9K w ■ {v -vk) - gci^rn^iv) ■ {v -VCs,) 



As explained in Sect. [2l the ML model tracks the membrane voltage v and a gating variable w. The 
constant Cm is the membrane capacitance, (p a timescale paramter, and I{t) an injected current. 

Spike generation depends crucially on the presence of voltage-sensitive ion channels that are 
permeable only to specific types of ions. The ML equations include just two ionic currents, here 
denoted calcium and potassium. The voltage response of ion channels is governed by the i(;-equation 
and the auxiliary functions Woo , , and rrioo , which have the form 



The function moo{v) models the action of the relatively fast calcium ion channels; vca. is the "rever- 
sal" (bias) potential for the calcium current and gca the corresponding conductance. The gating vari- 
able w and the functions Tu,{v) and Wao{v) model the dynamics of slower-acting potassium channels, 
with its own reversal potential vk and conductance gx ■ The constants viea.k and gieak characterize the 
"leakage" current that is present even when the neuron is in a "quiescent" state. The forms of rrioo , 
Tw , and Woo (as well as the values of the Vi) can be obtained by fitting data, or reduction from more 
biophysically-faithful models of Hodgkin-Huxley type (see, e.g., [9|). 

The precise parameter values used in this paper are summarized in Table [T] These are obtained 
from in. 
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